Analysis GEUVADIS LLS results

mydir.scripts <- "./scripts"
mydir.output <-  "./output"
mydir.data <- "./data"
mydir.output.perGene <- mydir.output
source(paste0(mydir.scripts, "/functions_wilcoxon_t_logistic_fdr.R"))

Results GEUVADIS and LLS

Data used

The LLS data was produced as part of the BBMRI project. Data can be requested for registered users via their website.

The GEUVADIS data has been pre-processed again in the same way as the LLS data, and analyses were rerun. The original GEUVADIS data can be obtained for example from the EBI.

The script used to install the spliceQTL package, obtain SNP and exon-level data, and produce objects needed, can be found on Armin Rauschenberger’s own [github](https://github.com/rauschenberger/spliceQTL/blob/master/vignettes/code.Rmd] page. Lines 1-131 refer to installing, preprocessing and running analyses. In the remainder results are explored. Some of the results in the article are included below.

Researchers willing to perform these analyses for other datasets are encouraged to use the chunks in that report.

Below we reproduce some of the graphs presented in our article, loading results obtained.

Results

The spliceQTL test was run in the same way as done for the analysis of the GEUVADIS data (with the original pre-processing). Here we load the results of that analysis.

vchrs <- 1:22
mtable <- matrix(0, nrow = length(vchrs), ncol = 4)
colnames(mtable) <- c("none", "onlyGEU", "onlyLLS", "both")
rownames(mtable) <- vchrs
mtable <- data.frame(mtable)
table.g <- table.l <- table.go <- table.lo <- NULL
mlist <- list()
for(chr in vchrs)
  {
   load(paste0(mydir.data, "/pvalNew.Geuvadis.chr", chr, ".RData")) # pvalue
   a <- pvalue; pvalue <- NA
   a <- data.frame(a, chr = rep(chr, length(a)))
   table.g <- rbind(table.g, a)
   load(paste0(mydir.data, "/pvalNew.LLS.chr", chr, ".RData")) # pvalue
   b <- pvalue; pvalue <- NA
   b <- data.frame(b, chr = rep(chr, length(b)))
   table.l <- rbind(table.l, b)

  # filter tests (trial)
  names <- intersect(rownames(a),rownames(b))
  a <- a[names, ]
  b <- b[names, ]

  # If desired, compute correlations between results of the 2 datasets
  if(FALSE){ 
  mtable$cor[ chr ] <- round(stats::cor(-log(a[, 1]), -log(b[, 1]), method="pearson"), digits=2)
  mtable$cor.test[ chr ] <- signif(stats::cor.test(-log(a[, 1]), -log(b[, 1]), method="pearson")$p.value, digits=2)
  }
  
  # significance
  a <- stats::p.adjust(a[, 1]) < 0.05
  b <- stats::p.adjust(b[, 1]) < 0.05
  mtable$none[ rownames(mtable) == chr ] <- sum(!a & !b, na.rm = TRUE)
  mtable$onlyGEU[ rownames(mtable) == chr ] <- sum(a & !b, na.rm = TRUE)
  mtable$onlyLLS[ rownames(mtable) == chr ] <- sum(!a & b, na.rm = TRUE)
  mtable$both[ rownames(mtable) == chr ] <- sum(a & b, na.rm = TRUE)
}
save(table.g, file=paste0(mydir.output, "/table_allpvalues_GEUV.RData"))
save(table.l, file=paste0(mydir.output, "/table_allpvalues_LLS.RData"))
# Tables with probes in the overlap are saved in the next chunk,
# after computing the FDR
save(mtable, file=paste0(mydir.output, "/table_cross_GEUV_LLS.RData"))

load(paste0(mydir.output, "/table_allpvalues_GEUV.RData")) # table.g
load(paste0(mydir.output, "/table_allpvalues_LLS.RData")) # table.l
load(paste0(mydir.output, "table_cross_GEUV_LLS.RData")) # mtable

mcut <- 0.05
dname <- "GEUVADIS"
mytable <- table.g
fdr <- NULL
for(chr in vchrs) 
{
  mname <- paste(dname, "Chr", chr)
  vfdr <- adj.bhfdr(mytable[mytable$chr == chr, 1], myvar = mname, cutoff = mcut)
  fdr <- c(fdr, vfdr)
}

mytable$fdr <- fdr
table.g <- mytable

###

dname <- "LLS"
mytable <- table.l
fdr <- NULL
for(chr in vchrs) 
{
  mname <- paste(dname, "Chr", chr)
  vfdr <- adj.bhfdr(mytable[mytable$chr == chr, 1], myvar = mname, cutoff = mcut)
  fdr <- c(fdr, vfdr)
}

mytable$fdr <- fdr
table.l <- mytable


names.o <- intersect(rownames(table.l),rownames(table.g))
table.lo <- table.l[names.o, ]
table.go <- table.g[names.o, ]
save(table.go, file=paste0(mydir.output, "/table_overlappvalues_GEUV.RData"))
save(table.lo, file=paste0(mydir.output, "/table_overlappvalues_LLS.RData"))

In the code chunk below plots of the FDRs per chromosome and dataset are produced. These make use of all results per dataset, not the results from the overlapping genes. Graphs are saved into a pdf file.

myn <- "GEUV"
myd <- table.g
pdf(paste0(mydir.output, "/plot_fdr_allChrs_", myn, ".pdf"), width = 8, height = 12)
par(mfrow = c(6, 4), mar = c(5, 4, 0.5, 0.5))
for(chr in vchrs)
{ 
   mf <- myd$fdr[myd$chr == chr]
   plot(sort(mf), type = "l", col = "blue", xlab = "genes", ylab = "FDR")
   mfs <- mf[mf <= 0.05]
   lines(sort(mfs), col = "red", lwd = 2)
   segments(0, 0.05, length(mf), 0.05, lty = "dashed", col = "gray")
   text(0, 0.7, labels = paste(length(mfs), "FDR<0.05"), pos = 4, cex = .8)
   text(0, 0.9, labels = paste("chr", chr), pos = 4, cex = .8)
}
dev.off()

myn <- "LLS"
myd <- table.l
pdf(paste0(mydir.output, "/plot_fdr_allChrs_", myn, ".pdf"), width = 8, height = 12)
par(mfrow = c(6, 4), mar = c(5, 4, 0.5, 0.5))
for(chr in vchrs)
{ 
   mf <- myd$fdr[myd$chr == chr]
   plot(sort(mf), type = "l", col = "blue", xlab = "genes", ylab = "FDR")
   mfs <- mf[mf <= 0.05]
   lines(sort(mfs), col = "red", lwd = 2)
   segments(0, 0.05, length(mf), 0.05, lty = "dashed", col = "gray")
   text(0, 0.7, labels = paste(length(mfs), "FDR<0.05"), pos = 4, cex = .8)
   text(0, 0.9, labels = paste("chr", chr), pos = 4, cex = .8)
}
dev.off()

Tables of significant results from the overlapping genes. First without multiple testing correction, and testing for association between the two test results. This is supplementary figure 10.

table.l <- table.l[rownames(table.l) %in% rownames(table.g), ]
table.g <- table.g[rownames(table.g) %in% rownames(table.l), ]
#all(rownames(table.g) == rownames(table.l))
dcols <- densCols(-log(table.g$a), -log(table.l$b))
plot(-log(table.g$a), -log(table.l$b), pch = 20, col = dcols,
     xlab = "Geuvadis", ylab = "LLS", main = "SpliceQTL p-values per gene")
abline(h = 3, lty = "dashed", col = "gray", lwd = 2)
abline(v = 3, lty = "dashed", col = "gray", lwd = 2)

Tables, no multiple testing correction. Different thresholds.

vth <- c(0.001, 0.005, 0.01)
for(th in vth)
{
   myt <- table(table.g$a <= th, table.l$b <= th)
   print(myt)
   print(chisq.test(myt)) # test of marginal independence, not quite what we want
   print(mcnemar.test(myt)) #test of symmetry
}

##        
##         FALSE  TRUE
##   FALSE 14915  1441
##   TRUE    421   603
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  myt
## X-squared = 2323.8, df = 1, p-value < 2.2e-16
## 
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  myt
## McNemar's chi-squared = 557.66, df = 1, p-value < 2.2e-16
## 
##        
##         FALSE  TRUE
##   FALSE 14168  1811
##   TRUE    590   811
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  myt
## X-squared = 2175.5, df = 1, p-value < 2.2e-16
## 
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  myt
## McNemar's chi-squared = 619.91, df = 1, p-value < 2.2e-16
## 
##        
##         FALSE  TRUE
##   FALSE 13213  2249
##   TRUE    830  1088
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  myt
## X-squared = 1954.2, df = 1, p-value < 2.2e-16
## 
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  myt
## McNemar's chi-squared = 653.04, df = 1, p-value < 2.2e-16

vth <- c(0.001, 0.005, 0.01)
for(th in vth)
{
   myt <- table(table.g$a <= th, table.l$b <= th)
   print(myt)
 #  print(chisq.test(myt)) # test of marginal independence, not quite what we want
 #  print(mcnemar.test(myt)) #test of symmetry
   myt <- 100*myt/sum(myt)
   print(chisq.test(myt)) # test of marginal independence, not quite what we want
   print(mcnemar.test(myt)) #test of symmetry
}

##        
##         FALSE  TRUE
##   FALSE 14915  1441
##   TRUE    421   603

## Warning in chisq.test(myt): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  myt
## X-squared = 9.0074, df = 1, p-value = 0.002689
## 
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  myt
## McNemar's chi-squared = 2.2127, df = 1, p-value = 0.1369
## 
##        
##         FALSE  TRUE
##   FALSE 14168  1811
##   TRUE    590   811

## Warning in chisq.test(myt): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  myt
## X-squared = 9.1674, df = 1, p-value = 0.002464
## 
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  myt
## McNemar's chi-squared = 2.628, df = 1, p-value = 0.105
## 
##        
##         FALSE  TRUE
##   FALSE 13213  2249
##   TRUE    830  1088

## Warning in chisq.test(myt): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  myt
## X-squared = 8.7048, df = 1, p-value = 0.003174
## 
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  myt
## McNemar's chi-squared = 2.8975, df = 1, p-value = 0.08872

The correlation between the two test results (on the -log scale) is 0.43081 for Pearson and 0.2708648 for Spearman.

Now with multiple testing correction.

table(table.go[, "fdr"] <= 0.05, table.lo[, "fdr"] <= 0.05)

##        
##         FALSE  TRUE
##   FALSE 13954  2137
##   TRUE    486   803

for(chr in vchrs) 
   {
     print(chr)
     print(table(table.go[table.go$chr == chr, "fdr"] <= 0.05, table.lo[table.lo$chr == chr, "fdr"] <= 0.05))
}

## [1] 1
##        
##         FALSE TRUE
##   FALSE  1526  194
##   TRUE     39   76
## [1] 2
##        
##         FALSE TRUE
##   FALSE   980  123
##   TRUE     40   61
## [1] 3
##        
##         FALSE TRUE
##   FALSE   773   88
##   TRUE     27   38
## [1] 4
##        
##         FALSE TRUE
##   FALSE   523   71
##   TRUE     22   32
## [1] 5
##        
##         FALSE TRUE
##   FALSE   661  102
##   TRUE     28   34
## [1] 6
##        
##         FALSE TRUE
##   FALSE   741   86
##   TRUE     33   37
## [1] 7
##        
##         FALSE TRUE
##   FALSE   669   93
##   TRUE     31   42
## [1] 8
##        
##         FALSE TRUE
##   FALSE   518   80
##   TRUE     16   28
## [1] 9
##        
##         FALSE TRUE
##   FALSE   564   74
##   TRUE     14   17
## [1] 10
##        
##         FALSE TRUE
##   FALSE   594   78
##   TRUE     13   26
## [1] 11
##        
##         FALSE TRUE
##   FALSE   879  135
##   TRUE     26   44
## [1] 12
##        
##         FALSE TRUE
##   FALSE   820   93
##   TRUE     36   49
## [1] 13
##        
##         FALSE TRUE
##   FALSE   265   22
##   TRUE      7   12
## [1] 14
##        
##         FALSE TRUE
##   FALSE   486   65
##   TRUE     13   28
## [1] 15
##        
##         FALSE TRUE
##   FALSE   403   83
##   TRUE     23   33
## [1] 16
##        
##         FALSE TRUE
##   FALSE   603  141
##   TRUE     14   36
## [1] 17
##        
##         FALSE TRUE
##   FALSE   801  191
##   TRUE     27   53
## [1] 18
##        
##         FALSE TRUE
##   FALSE   222   27
##   TRUE     10   14
## [1] 19
##        
##         FALSE TRUE
##   FALSE  1050  226
##   TRUE     39   81
## [1] 20
##        
##         FALSE TRUE
##   FALSE   436   62
##   TRUE     11   15
## [1] 21
##        
##         FALSE TRUE
##   FALSE   134   31
##   TRUE      4   11
## [1] 22
##        
##         FALSE TRUE
##   FALSE   306   72
##   TRUE     13   36

Now Manhattan plots are produced, saved into pdfs for 5 consecutive chromosomes at a time. These correspond to supplementary figures 11–15.

# This uses the overlap
dname <- "GEUVADIS and LLS"
pdf(paste0(mydir.output, "/plot_manhattan_both_1.pdf"), width = 8, height = 12)
par(mfrow = c(5, 1), mar = c(2, 3, 3, 1))
for(chr in vchrs[1:5]) 
{
   mname <- paste(dname, "Chr", chr)
   a <- table.go[table.go$chr == chr, "fdr"]
   b <- table.lo[table.lo$chr == chr, "fdr"]
   graphics::plot.new()
   graphics::plot.window(xlim=c(0.5, length(a)+0.5), ylim=c(-1,1)*max(-log10(c(a, b))))
   graphics::box()
   graphics::axis(side=1)
   graphics::axis(side=2)
   graphics::segments(x0=seq_along(a), y0=0, y1=-log10(a), col="darkred")
   graphics::segments(x0=seq_along(b), y0=0, y1=log10(b), col="darkblue")
   graphics::abline(h=0, col="grey")
   title(paste("Chr", chr, "FDR - Geuvadis (red) and LLS (blue)"))
}
dev.off()

pdf(paste0(mydir.output, "/plot_manhattan_both_2.pdf"), width = 8, height = 12)
par(mfrow = c(5, 1), mar = c(2, 3, 3, 1))
for(chr in vchrs[6:10]) 
{
   mname <- paste(dname, "Chr", chr)
   a <- table.go[table.go$chr == chr, "fdr"]
   b <- table.lo[table.lo$chr == chr, "fdr"]
   graphics::plot.new()
   graphics::plot.window(xlim=c(0.5, length(a)+0.5), ylim=c(-1,1)*max(-log10(c(a, b))))
   graphics::box()
   graphics::axis(side=1)
   graphics::axis(side=2)
   graphics::segments(x0=seq_along(a), y0=0, y1=-log10(a), col="darkred")
   graphics::segments(x0=seq_along(b), y0=0, y1=log10(b), col="darkblue")
   graphics::abline(h=0, col="grey")
   title(paste("Chr", chr, "FDR - Geuvadis (red) and LLS (blue)"))
}
dev.off()

pdf(paste0(mydir.output, "/plot_manhattan_both_3.pdf"), width = 8, height = 12)
par(mfrow = c(5, 1), mar = c(2, 3, 3, 1))
for(chr in vchrs[11:15]) 
{
   mname <- paste(dname, "Chr", chr)
   a <- table.go[table.go$chr == chr, "fdr"]
   b <- table.lo[table.lo$chr == chr, "fdr"]
   graphics::plot.new()
   graphics::plot.window(xlim=c(0.5, length(a)+0.5), ylim=c(-1,1)*max(-log10(c(a, b))))
   graphics::box()
   graphics::axis(side=1)
   graphics::axis(side=2)
   graphics::segments(x0=seq_along(a), y0=0, y1=-log10(a), col="darkred")
   graphics::segments(x0=seq_along(b), y0=0, y1=log10(b), col="darkblue")
   graphics::abline(h=0, col="grey")
   title(paste("Chr", chr, "FDR - Geuvadis (red) and LLS (blue)"))
}
dev.off()

pdf(paste0(mydir.output, "/plot_manhattan_both_4.pdf"), width = 8, height = 12)
par(mfrow = c(5, 1), mar = c(2, 3, 3, 1))
for(chr in vchrs[16:20]) 
{
   mname <- paste(dname, "Chr", chr)
   a <- table.go[table.go$chr == chr, "fdr"]
   b <- table.lo[table.lo$chr == chr, "fdr"]
   graphics::plot.new()
   graphics::plot.window(xlim=c(0.5, length(a)+0.5), ylim=c(-1,1)*max(-log10(c(a, b))))
   graphics::box()
   graphics::axis(side=1)
   graphics::axis(side=2)
   graphics::segments(x0=seq_along(a), y0=0, y1=-log10(a), col="darkred")
   graphics::segments(x0=seq_along(b), y0=0, y1=log10(b), col="darkblue")
   graphics::abline(h=0, col="grey")
   title(paste("Chr", chr, "FDR - Geuvadis (red) and LLS (blue)"))
}
dev.off()

pdf(paste0(mydir.output, "/plot_manhattan_both_5.pdf"), width = 8, height = 12)
par(mfrow = c(5, 1), mar = c(2, 3, 3, 1))
for(chr in vchrs[21:22]) 
{
   mname <- paste(dname, "Chr", chr)
   a <- table.go[table.go$chr == chr, "fdr"]
   b <- table.lo[table.lo$chr == chr, "fdr"]
   graphics::plot.new()
   graphics::plot.window(xlim=c(0.5, length(a)+0.5), ylim=c(-1,1)*max(-log10(c(a, b))))
   graphics::box()
   graphics::axis(side=1)
   graphics::axis(side=2)
   graphics::segments(x0=seq_along(a), y0=0, y1=-log10(a), col="darkred")
   graphics::segments(x0=seq_along(b), y0=0, y1=log10(b), col="darkblue")
   graphics::abline(h=0, col="grey")
   title(paste("Chr", chr, "FDR - Geuvadis (red) and LLS (blue)"))
}
dev.off()

Scatterplots of raw p-values, as well as FDR-corrected p-values, per chromosome are produced below. Results are saved into pdf files.

# pvalues
dname <- "GEUVADIS and LLS"
pdf(paste0(mydir.output, "/plot_pvalues_GEUV_LLS_perChr.pdf"), width = 10, height = 5)
par(mfrow = c(1, 2))
for(chr in vchrs) 
{
   mname <- paste(dname, "Chr", chr)
   a <- -log10(table.go[table.go$chr == chr, 1])
   b <- -log10(table.lo[table.lo$chr == chr, 1])
   mylims <- range(c(a, b))
   dcols <- densCols(a, b)
   plot(a, b, col = dcols, pch = 20, main = mname, xlab = "Geuvadis -log10 pvalue", ylab = "LLS -log10 pvalue")
   segments(0, 0, mylims[2], mylims[2], lty = "dashed", col = "gray", lwd = 2)
   segments(0, -log10(0.05), mylims[2], -log10(0.05), lty = "dashed", col = "gray", lwd = 2)
   segments(-log10(0.05), 0, -log10(0.05), mylims[2], lty = "dashed", col = "gray", lwd = 2)
}
dev.off()

# fdr
pdf(paste0(mydir.output, "/plot_fdrs_GEUV_LLS_perChr.pdf"), width = 10, height = 5)
for(chr in vchrs) 
{
   mname <- paste(dname, "Chr", chr)
   a <- -log10(table.go[table.go$chr == chr, "fdr"])
   b <- -log10(table.lo[table.lo$chr == chr, "fdr"])
   mylims <- range(c(a, b))
   dcols <- densCols(a, b)
   plot(a, b, col = dcols, pch = 20, main = mname, xlab = "Geuvadis -log10 FDR", ylab = "LLS -log10 FDR")
   segments(0, 0, mylims[2], mylims[2], lty = "dashed", col = "gray", lwd = 2)
   segments(0, -log10(0.05), mylims[2], -log10(0.05), lty = "dashed", col = "gray", lwd = 2)
   segments(-log10(0.05), 0, -log10(0.05), mylims[2], lty = "dashed", col = "gray", lwd = 2)
}
dev.off()

Scatterplots of the FDR-corrected p-values of both datasets are also produced, and saved into a pdf file.

# fdr
pdf(paste0(mydir.output, "/plot_fdrs_bothDatasets.pdf"), width=8, height = 8)
   a <- -log10(table.go[, "fdr"])
   b <- -log10(table.lo[, "fdr"])
   mylims <- range(c(a, b))
   dcols <- densCols(a, b)
   plot(a, b, col = dcols, pch = 20, main = "All overlapping genes", xlab = "Geuvadis -log10 FDR", ylab = "LLS -log10 FDR")
   segments(0, 0, mylims[2], mylims[2], lty = "dashed", col = "gray", lwd = 2)
   segments(0, -log10(0.05), mylims[2], -log10(0.05), lty = "dashed", col = "gray", lwd = 2)
   segments(-log10(0.05), 0, -log10(0.05), mylims[2], lty = "dashed", col = "gray", lwd = 2)
dev.off()

Below scatterplots of the raw and FDR-corrected p-values of both datasets are produced, as well as of the p-values only. Both are saved as a pdf file.

pdf(paste0(mydir.output, "/plot_pvals_fdr_Geuv_LLS.pdf"), width = 5, height = 10)
par(mfrow = c(2, 1))
dcols <- densCols(-log10(table.g$a), -log10(table.l$b))
plot(-log10(table.g$a), -log10(table.l$b), pch = 20, col = dcols, cex = 0.7,
     xlab = "Geuvadis", ylab = "LLS", main = "P-value per gene")
   mylims <- range(c(-log10(table.g$a), -log10(table.l$b)))
   segments(0, -log10(0.001), mylims[2], -log10(0.001), lty = "dashed", col = "gray", lwd = 2)
   segments(-log10(0.001), 0, -log10(0.001), mylims[2], lty = "dashed", col = "gray", lwd = 2)
   a <- -log10(table.go[, "fdr"])
   b <- -log10(table.lo[, "fdr"])
   mylims <- range(c(a, b))
   dcols <- densCols(a, b)
   plot(a, b, col = dcols, pch = 20, main = "FDR per gene", cex = 0.7, xlab = "Geuvadis", ylab = "LLS")
   segments(0, 0, mylims[2], mylims[2], lty = "dashed", col = "gray", lwd = 2)
   segments(0, -log10(0.05), mylims[2], -log10(0.05), lty = "dashed", col = "gray", lwd = 2)
   segments(-log10(0.05), 0, -log10(0.05), mylims[2], lty = "dashed", col = "gray", lwd = 2)

dev.off()

pdf(paste0(mydir.output, "/plot_pvals_only_Geuv_LLS.pdf"), width = 6, height = 6)
dcols <- densCols(-log10(table.g$a), -log10(table.l$b))
plot(-log10(table.g$a), -log10(table.l$b), pch = 20, col = dcols, cex = 0.7,
     xlab = "Geuvadis", ylab = "LLS", main = "SpliceQTL results per gene")
   mylims <- range(c(-log10(table.g$a), -log10(table.l$b)))
   segments(0, -log10(0.001), mylims[2], -log10(0.001), lty = "dashed", col = "gray", lwd = 2)
   segments(-log10(0.001), 0, -log10(0.001), mylims[2], lty = "dashed", col = "gray", lwd = 2)
dev.off()

Make selection

Now we select genes that are statistically significant in either dataset, or both.

mcut <- 0.01
table.go$sel <- table.go[, "fdr"] <= mcut
table.lo$sel <- table.lo[, "fdr"] <= mcut
table(table.go$sel, table.go$chr)

##        
##            1    2    3    4    5    6    7    8    9   10   11   12   13   14
##   FALSE 1738 1147  875  608  785  844  786  617  647  681 1034  932  291  561
##   TRUE    97   57   51   40   40   53   49   25   22   30   50   66   15   31
##        
##           15   16   17   18   19   20   21   22
##   FALSE  504  764 1013  261 1305  503  171  389
##   TRUE    38   30   59   12   91   21    9   38

table(table.lo$sel, table.lo$chr)

##        
##            1    2    3    4    5    6    7    8    9   10   11   12   13   14
##   FALSE 1629 1078  843  579  725  801  721  566  610  640  943  885  278  517
##   TRUE   206  126   83   69  100   96  114   76   59   71  141  113   28   75
##        
##           15   16   17   18   19   20   21   22
##   FALSE  466  680  923  245 1189  466  154  362
##   TRUE    76  114  149   28  207   58   26   65

#table.go.sel$sel <- table.go$fdr <= mcut
#table.lo.sel$sel <- table.lo$fdr <= mcut
for(xch in vchrs)
{ 
 print(paste("Chr", xch))
   print(table(table.go[table.go$chr == xch, "sel"], table.lo[table.lo$chr == xch, "sel"]))
}

## [1] "Chr 1"
##        
##         FALSE TRUE
##   FALSE  1595  143
##   TRUE     34   63
## [1] "Chr 2"
##        
##         FALSE TRUE
##   FALSE  1058   89
##   TRUE     20   37
## [1] "Chr 3"
##        
##         FALSE TRUE
##   FALSE   814   61
##   TRUE     29   22
## [1] "Chr 4"
##        
##         FALSE TRUE
##   FALSE   564   44
##   TRUE     15   25
## [1] "Chr 5"
##        
##         FALSE TRUE
##   FALSE   707   78
##   TRUE     18   22
## [1] "Chr 6"
##        
##         FALSE TRUE
##   FALSE   775   69
##   TRUE     26   27
## [1] "Chr 7"
##        
##         FALSE TRUE
##   FALSE   704   82
##   TRUE     17   32
## [1] "Chr 8"
##        
##         FALSE TRUE
##   FALSE   554   63
##   TRUE     12   13
## [1] "Chr 9"
##        
##         FALSE TRUE
##   FALSE   601   46
##   TRUE      9   13
## [1] "Chr 10"
##        
##         FALSE TRUE
##   FALSE   630   51
##   TRUE     10   20
## [1] "Chr 11"
##        
##         FALSE TRUE
##   FALSE   928  106
##   TRUE     15   35
## [1] "Chr 12"
##        
##         FALSE TRUE
##   FALSE   860   72
##   TRUE     25   41
## [1] "Chr 13"
##        
##         FALSE TRUE
##   FALSE   272   19
##   TRUE      6    9
## [1] "Chr 14"
##        
##         FALSE TRUE
##   FALSE   508   53
##   TRUE      9   22
## [1] "Chr 15"
##        
##         FALSE TRUE
##   FALSE   449   55
##   TRUE     17   21
## [1] "Chr 16"
##        
##         FALSE TRUE
##   FALSE   668   96
##   TRUE     12   18
## [1] "Chr 17"
##        
##         FALSE TRUE
##   FALSE   900  113
##   TRUE     23   36
## [1] "Chr 18"
##        
##         FALSE TRUE
##   FALSE   240   21
##   TRUE      5    7
## [1] "Chr 19"
##        
##         FALSE TRUE
##   FALSE  1155  150
##   TRUE     34   57
## [1] "Chr 20"
##        
##         FALSE TRUE
##   FALSE   458   45
##   TRUE      8   13
## [1] "Chr 21"
##        
##         FALSE TRUE
##   FALSE   151   20
##   TRUE      3    6
## [1] "Chr 22"
##        
##         FALSE TRUE
##   FALSE   350   39
##   TRUE     12   26

# Here select ids to make plots in the server
glist <- tbsel <- NULL
mlo <- 10^(-2)
mlo.l <- 10^(-2.5)
mup <- 10^(-1)
for(xch in 1:22)
{ 
   tbg <- table.go[(table.go$chr == xch), ]
   colnames(tbg)[3] <- "fdr.g"
   tbl <- table.lo[(table.lo$chr == xch), ]
   colnames(tbl)[3] <- "fdr.l"
   tbl$chr <- NULL
   tb2 <- cbind(tbg, tbl)
   tbsel1 <- tb2[(tb2$fdr.g <= mlo)&(tb2$fdr.l > mup), ]
   g1 <- rownames(tb2[(tb2$fdr.g <= mlo)&(tb2$fdr.l > mup), ] )
   tbsel1 <- rbind(tbsel1, tb2[(tb2$fdr.l <= mlo.l)&(tb2$fdr.g > mup), ])
   g2 <- rownames(tb2[(tb2$fdr.l <= mlo.l)&(tb2$fdr.g > mup), ] )
   tbsel1 <- rbind(tbsel1, tb2[(tb2$fdr.g <= mlo)&(tb2$fdr.l <= mlo), ])
   g3 <- rownames(tb2[(tb2$fdr.g <= mlo)&(tb2$fdr.l <= mlo), ] )
   vg <- c(rep("g", length(g1)), rep("l", length(g2)), rep("both", length(g3)))
   glist <- rbind(glist, matrix(c(g1, g2, g3, vg), nrow = length(vg), ncol = 2))
   tbsel1 <- cbind(tbsel1, vg)
   tbsel <- rbind(tbsel, tbsel1)
}
save(tbsel, file = paste0(mydir.output, "/table_sel_geuv_lls.RData")) # tbsel

Table of selections per chr arm.

table(tbsel$chr, tbsel$vg)

##     
##      both   g   l
##   1    63  24 101
##   2    37  10  50
##   3    22  15  37
##   4    25   8  30
##   5    22  13  50
##   6    27  20  44
##   7    32  15  45
##   8    13   4  39
##   9    13   6  34
##   10   20   6  34
##   11   35  10  62
##   12   41  18  54
##   13    9   4  14
##   14   22   5  32
##   15   21  10  40
##   16   18   6  73
##   17   36  14  76
##   18    7   5   9
##   19   57  19 111
##   20   13   5  37
##   21    6   2  11
##   22   26   8  27

Technical information

sessionInfo()

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=nl_NL.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=nl_NL.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=nl_NL.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=nl_NL.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.29      R6_2.5.0           jsonlite_1.7.2     magrittr_2.0.1    
##  [5] evaluate_0.14      KernSmooth_2.23-17 highr_0.9          stringi_1.6.2     
##  [9] rlang_0.4.12       jquerylib_0.1.4    bslib_0.2.5.1      rmarkdown_2.9     
## [13] tools_3.6.3        stringr_1.4.0      xfun_0.24          yaml_2.2.1        
## [17] fastmap_1.1.0      compiler_3.6.3     htmltools_0.5.2    knitr_1.33        
## [21] sass_0.4.0